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ABSTRACT 



The multi-spin coding of the Monte Carlo simulation of the three-state Potts model 
d ' on the simple cubic lattice is presented. The ferromagnetic (F) model, the antiferro- 

magnetic (AF) model, and the random mixture of the F and AF couplings are treated. 
The multi-spin coding technique is also applied to the block-spin transformation. The 
block-spin transformation of the F Potts model is simply realized by the majority rule, 
whereas the AF three-state Potts model is transformed to the block spin having a 



, six-fold symmetry. 



1. Introduction 

a' 



Monte Carlo simulations are used as standard techniques to investigate statistical mechan- 
ical properties of many-body systems. In treating large systems, especially near the critical 
points and at very low temperatures, we often encounter slow dynamics. To overcome such 
slow dynamics, the development of fast algorithms is demanded. To gather more infor- 
mation from a single simulation is another direction of effort for algorithmic improvement. 
The cluster-flip Monte Carlo methodJlH the histogram method,! and the multicanonical 
simulationcl are examples of the recent progress. 

For the simulations of the Ising model, where only one bit is required for storing the 
information of a single spin, a computer word can store several spins. Based on this fact, the 
multi-spin coding technique has been successfully developedB&Qfl Among them, the idea 
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of the coding by Bhanot et a/.i and its variations are especially useful.0'00Jll In these 
methods, just a single bit is assigned to each spin. Instead of storing spins at different sites 
in the same lattice in a word, as in the case of former realizations of the multi spin coding, 
Ising spins at the same site of several independent systems (of the same lattice structure) are 
stored in a word. The spins for 32 (64) systems are updated simultaneously in the case of 32- 
(64-) bit machine, with a single random-number sequence; for that to be achieved, all the 
operations are executed by logical commands. As a result, the computation time is reduced 
remarkably. One can simulate either systems with different parameters, for example, the 
temperature, the external field, etc., or systems under the same parameters with different 
random number sequencea^-EYeri for the latter purpose, one needs to generate only one 
random-number sequence In the case of the simulations of random systems, one 
can simulate systems with different configurations in parallel.0 Since an average over a 
large number of random configurations should be taken, the multi-spin coding technique is 
particularly effective for random systems. We should note that the multi-spin coding has 
also been applied to the Monte Carlo simulations of the Ising model on quasicrystalsOEl 

The Monte Carlo renormalization group (MCRG) is a powerful tool to analyze criti- 
cal phenomena.0 In MCRG method, one performs the block-spin transformations at each 
Monte Carlo sweep. The block-spin transformation makes a block spin from b D original 
spins, with the scale factor b. The present authors have successfully used the multi-spin 
coding to the block-spin transformation of the three-dimensional (3D) Ising model.t!-!ffi3 

The idea of the multi-spin coding can also be applied to other models with discrete 
symmetries. The g-state Potts model is one example. The Potts model has various in- 
teresting properties, depends The order of the phase transitions of the ferromagnetic (F) 
Potts model depend on the dimensionality, D, and the number of the state, q. In 3D, the 
F Potts model shows a first-order transition, whereas the antiferromagnetic (AF) model 
shows a second-order transition. The random mixture of F and AF couplings gives rise to 
new problems: The rounding of first-order transition in weak ferromagnetic region is one 
of interesting problems; one may also ask the universality of critical exponents in this case. 
The Potts glass phase will be another problem of interest. The zero-temperature transition 
of the three-state Potts glass has been suggested in 3D. 

In this paper, we present the multi-spin coding of the three-state Potts model on the 
simple cubic lattice. For the coupling, we consider both cases of the F and AF couplings. 
The random mixture of the F and AF couplings are also treated. We also present the 
realization of the block-spin transformation using multi-spin coding technique. The appro- 
priate choice of the block-spin transformation is essential to extract the critical properties of 
the order parameter. We employ quite different block-spin transformations for the F Potts 
model and the AF Potts model: The block-spin transformation of the F model is simply 
realized by the majority rule, whereas the block spin of the AF three-state Potts model has 
the same symmetry as the six-clock model. 

2. Multi-Spin Coding 

We are concerned with the three-state Potts Model, whose Hamiltonian is given as 

w = - E J iMst , (i) 

<i,3> 
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where each spin Si can take three states, 0, 1 and 2. Here a variable Jij may take either 
+ J or — J generally. In the following, we consider only the nearest-neighbor interactions. 

We employ the multi-spin coding algorithmpliHI which was originally used for the Ising 
model. In contrast to the Ising model, two bits are required to represent a Potts spin; 
the three states, 0, 1 and 2, are represented by (00), (01), and (10) in binary-number 
representation. Therefore, we assign one word for the upper bit and another word for the 
lower bit, as is seen in the following. It should be noted that the calculation of the order 
parameter becomes quite simple in this two-word representation, because the number of 1 
(2) spin is counted just by summing up the lower- (upper-) bit throughout the lattice. 

To update a spin in one of the three states, one has to choose either of the rest two 
states as a trial state. This procedure is executed by the following code: 
iscoin=ira(la) 

isnewO=iand(not (isoldO) ,iscoin) 
isnewl=not (ior (isoldl , iscoin) ) 

The old spin is represented by (isoldl, isoldO), with the former word representing the 
upper bit and the latter the lower bit, respectively. Here, ira is a random number, where 
each bit takes or 1 with the probability of 1/2. Then, the trial state, also represented 
by two words, (isnewl, isnewO), is selected as one of the other two states with an equal 
probability for every case of the old states (00), (01), and (10) . 

In the Metropolis algorithm, one needs to count the local energy change due to the 
above trial flip of a picked single spin. The following code is for calculating the energy 
change (in unit of J) at a bond: 

neqold = ior(ieor(isold0,is0) ,ieor(isoldl,isl)) 
neqnew = ior (ieor (isnewO, isO) ,ieor (isnewl, isl)) 
ideO = iand (neqold, neqnew) 
idel = not (neqnew) 

Two words (isl, isO) represent the nearest-neighbor spin in concern. The energy change 
is also given by the binary-number representation in two words, idel and ideO. The above 
code is for F Potts model. For AF model, the last line should be modified as follows: 

idel = not (neqold) 
And for the random-bond model, the same line should be read as 

idel = ior (iand (jbond, not (neqnew) ), not (ior (jbond, neqold))) 

where jbond should either be 1 for the ferromagnetic bond or for the antiferromagnetic 
bond. Since the coordination number is six for the simple cubic lattice, the local energy 
change W takes an integer value between —6 and 6 in units of J. Thus, four words are 
necessary to represent W in the multi-spin coding. It should be compared to the Ising case, 
where three words are enough to represent the local energy change. 

In using the random-number sequence, we apply the trick proposed by MichaelJIS rather 
than the original one by Bhanot et aZ.i That is, we use the precalculated table for the 
transition probability, min[l, e~@ JW ]; in the original algorithm, on the other hand, the 
transition probability table is updated dynamically. In case of using the precalculated 
table, a large table is required to ensure the accuracy in temperature; one can, however, 
use a very large table easily in computers these days, so that it is not considered as a 
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drawback. On the other hand, using a small table which is dynamically updated may cause 
an uncontrolable temperature fluctuation. 

We use a dummy variable X, which takes the integer value from to 6. We prepare the 
table so that the probability to get the value of X is given by 
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where K = (5 J and r is a random number which uniformly takes the value between and 
1. Then, the rule of the spin update will be as follows: If (-W) + X + 8 is greater than 8, 
then the spin will be flipped. 

It should be emphasized that the shuffling of the precalculated A-table is effective for 
simulating different samples with the same parameters by a single random-number sequence. 
The statistical dependence of random sequences in conjunction with the shuffling of the table 
has been recently discussed .0 

The FORTRAN code for the Metropolis method is given in Appendix, which is for the 
F Potts model. In our multi-spin coding, X is represented by three words, (jxtb2, jxtbl, 
jxtbO). We use the integer random number, irb, which takes the value between and 
nword-1. The decrease in the local energy at the trial flip, (— W) + 8, is calculated to be 
(isum3, isum2, isuml, isumO) in binary-number representation. 

The spin update is executed by the logical commands; if lip is the flag for the flip. The 
updated spin and the resultant energy difference are given by ispO and iengd, respectively. 
With slight modifications, which are commented out with we can also treat the AF 
model and the ±J model; the information of the configuration of ±J bonds should be 
included in jbond, which in general depends on the lattice point la. For the convenience 
of vectorization, we divide the lattice into two interpenetrating sublattices. 

The linear size of the system is given by nx; the number of the spins on each sublattice is 
given by nla2=nx**3/2. The periodic boundary conditions are employed, and informations 
on the neighboring lattice points are provided in the tables, jx, jyr, jyl, jzr and jzl. 

3. Block-Spin Transformation 

3.1. Ferromagnetic Order Parameter 

For the three-state F Potts model, it is convenient to express the order parameter as a two- 
dimensional vector. Three states 0, 1, and 2 are expressed by the vectors directing ±2-7r/3 
from each other as shown in the inset of Fig. 1. Let us consider the case of the block-spin 
transformation from eight spins; that is, the scale factor b is 2. Possible vectors obtained 
by the vector sum of eight Potts spins are shown in Fig. 1. From each vector, we make a 
block spin which also takes three states as does the original Potts spin; we choose the Potts 
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Fig. 1. Possible vectors obtained by the vector sum of eight F Potts spins. 

state directing closest to the vector as the block-spin state. The rule for determination of 
the block spin is illustrated in Fig. 1, where the separation is shown with solid lines. At the 
border, the random number will be used to determine which state is to be chosen. 

One would find that the block-spin transformation described above is, in fact, equivalent 
to the simple majority rule; in other words, the three-state Potts spin is mapped onto the 
three-state Potts spin by the majority rule transformation. 

The essential part of the FORTRAN code for the block-spin transformation is given as 
follows: 

c 

c block-spin transformation 

c ferro 3-state Potts to 3-state Potts 

c 

non02x = ior(isumxO, isumxl) 
nsumx2 = not(isumx2) 
non02y = ior(isumyO, isumyl) 
nsumy2 = not(isumy2) 

jO = iand(iand(isumxO , isumxl) , isumyl) 
jl = iand(iand(isumyO , isumyl) , isumxl) 

c 

iran = ir(la) 

iranx = iand(isumx2 , iran) 

irany = iand(isumy2 , not (iran) ) 

c 

isO = iand(isumx2,iand(nsumy2,non02y)) 

isO = ior(is0,isumx3) 

isO = ior(is0,iand(isumx2,non02x)) 
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isO = ior(is0,iand(isumx2,irany)) 
isO = ior(is0,iand(nsumy2,iranx)) 
ibnewO(la) = ior (isO , iand( jO , iran) ) 

c 

isl = iand(isumy2,iand(nsumx2,non02x)) 

isl = ior (isl , isumy3) 

isl = ior (isl , iand(isumy2 ,non02y) ) 

isl = ior (isl , iand(isumy2 , iranx) ) 

isl = ior(isl,iand(nsumx2,irany)) 

ibnewl(la) = ior (isl , iand(j 1 ,not (iran) ) ) 

Here, the sum of 0-bit of eight Potts spins is expressed by (isumx3, isumx2, isumxl, isumxO). 
Similarly, the sum of 1-bit is given by (isumy3, isumy2, isumyl, isumyO). Two words, 
(ibnewl, ibnewO), gives the block spin. The random number ir is used for tie breaking at 
the border. To make the MCRG analysis of the F order parameter, we repeat the same 
block-spin transformation step by step. 



3.2. Antiferromagnetic Order Parameter 

For the AF order parameter, one should take account of the two-sublattice structure. Three 
components of the order parameter are defined as follows: 

^(E^-E^)/^, (2) 

ieA jeB 

where a denotes one of the three states 0, 1 and 2, and A and B stand for two sublattices. 
Three states 0, 1, and 2 for the sublattices A and B can be expressed by two-dimensional 
vectors as shown in the inset of Fig. 2. It should be noted here that each pair of the 
neighboring vectors make an angle of 7r/3; as a result, the order-parameter space in this 
case has a six-fold symmetry in contrast to that of the F Potts model. This symmetry, 
Z(6), is the same as that of the six-clock model. Based on this observation, we make block 
spins so as to preserve this six-fold symmetry; more specifically, as the block spin we choose 
one of the possible states of the six-clock spin which directs closest to the vector sum of 
the Potts spins in a block. Possible vectors obtained by the vector sum of eight Potts spins 
are illustrated in Fig. 2. At the border, a block spin should be chosen in equal probability 
by the random number. At the center, on the other hand, a block spin is chosen with the 
probability of 1/6. 

We have chosen the primary direction of the six-clock spin as shown in Fig. 2. The 

ered state of the AF three-state Potts model 
Namely, one sublattice is occupied by one of 
the three states, while the other sublattice is randomly occupied by the remaining two 
states. The present choice of the spin direction takes this type of order into account. 

The block-spin transformation described above for the AF order cannot be expressed 
as a simple majority rule in contrast to the F case. Therefore, we must perform the vector 
summation directly. The multi-spin coding technique is applicable even in this case, since 
all the operations can be executed only with integer variables space, when the oblique 
coordinates are used for expressing the order-parameter space. The six states, 0, 1, 2, 3, 



reason is as follows: We expect that theo 
has a broken-sublattice-symmetry order.BH£2 
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Fig. 2. Possible vectors obtained by the vector sum of eight AF Potts spins. 



4, 5 in Fig. 2, are expressed by (000), (001), (101), (100), (111), (011) respectively in our 
binary-number representation. We use the basis vectors e x and e y shown in Fig. 2. Then, 
0, 1 and 2 states on the sublattice A are expressed by r^4, e x + r^4 and e y + r^4, respectively, 
where = (e x + e y )/3. On the other hand, 0, 1 and 2 states on the sublattice B are 
expressed by e x + e y + re, e y + r# and e x + re, respectively, where rs = 2(e x + e y )/3. 

The FORTRAN code of the block-spin transformation from the AF Potts spin to the 
six-clock spin is given as follows: 

c 

c block spin transformation 

c antiferro 3-state Potts to 6-clock 

c 

c x>0 (icx = 0) or x=0 (icx = 1) or x<0 (icx = 2) 

c isumx >=< 4 

c 

iwk = ior (isumxO , isumx 1) 
icxO = iand(isumx2 ,not (iwk) ) 
icxl = not(ior(isumx2,isumx3)) 

c 

c y>0 (icy = 2) or y=0 (icy = 1) or y<0 (icy = 0) 

c isumy >=< 4 

c 

iwk = ior (isumyO , isumyl) 

icyO = iand(isumy2,not (iwk) ) 

icyl = iand(ior (isumy2, isumy3) ,not(icyO)) 

c 

c isumy + isumx (0-16) 
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icrO = iand(isumyO, isumxO) 
isumCO = ieor (isumyO, isumxO) 
iwk = ieor (isumyl , isumxl) 

icrl = ieor (iand(isumyl , isumxl) ,iand(icrO, iwk)) 
isumCl = ieor (icrO, iwk) 
iwk = ieor(isumy2,isumx2) 

icr2 = ieor(iand(isumy2,isumx2) , iandQcrl , iwk) ) 
isumC2 = ieor (icrl , iwk) 
iwk = ieor(isumy3,isumx3) 
isumC3 = ieor (icr2, iwk) 

isumC4 = ieor (iand(isumy3, isumx3) , iand(icr2 , iwk) ) 

c 

c y+x>0 (icm = 2) or y+x=0 (icm = 1) or y+x<0 (icm = 0) 
c isumy + isumx >=< 8 



c 



c 



iwk = ior(isumC0,ior(isumCl,isumC2)) 

icmO = iand(isumC3,not (iwk) ) 

icml = iand(ior (isumC3, isumC4) ,not(icmO)) 



c 

c if border 



iran = ira(la) 

icb2 = ior (icml , iand(icm0, iran) ) 
iwkl = ior(icxl,iand(icxO,iran)) 
iwk2 = ior (icy 1 , iand(icy0, iran) ) 
icbl = iand(iwkl , iwk2) 
icbO = not (ior (iwkl , iwk2) ) 



c 

c if center 
c 



icent = iand(icx0 , iand(icy0 , icmO) ) 
nocent = not (icent) 

c 

ic2 = irb(la) 

icl = not(ior(irtl(la) ,iran)) 
irtl(la) = icl 

icO = iand(not(irtO(la)) ,iran) 
irtO(la) = icO 

c 

isblk2(la) = ior (iand(nocent , icb2) , iand (icent , ic2) ) 
isblkl(la) = ior (iand(nocent , icbl) ,iand(icent, icl)) 
isblkO(la) = ior (iand(nocent , icbO) , iand(icent , icO) ) 

The sum of the x- and y-components of eight Potts spins are given by isumx and isumy. 
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The conditions needed for the determination of the block spin are represented in terms of 
the conditions on isumx, isumy and isumx+isumy, which are illustrated in Fig. 3. According 



Fig. 3. The conditions needed for the determination of block spin are illustrated. 

to these conditions, the three words to denote the six-clock spin, (isblk2, isblkl, isblkO), 
are determined by logical commands. The choice of the state at the border is simply made 
by the random number ira with the probability of 1/2. At the center, on the other hand, 
we use two random numbers, ira and irb: First, isblk2 is determined by the random 
number irb; next, the same procedure is employed for the determination of the three states 
represented by isblkl and isblkO as we used in choosing a trial state in the Monte Carlo 
update. The information of the selected state, (irtl, irtO), is kept up to the next time. 
With this procedure, we can pick up one state from six with the equal probabilities. 

Once the block spin having a six-fold symmetry is obtained, the next step will be 
the transformation from the six-clock spin to the block spin of the same symmetry. The 
illustration of this transformation is given in Fig. 4. Ties at the border and the center are 
treated similarly as before. We can apply the multi-spin coding technique also to the block- 
spin transformation from the six-clock spin to the six-clock block spin as before. Although 
we do not show the FORTRAN code for this transformation, it is worthwhile to make one 
comment: This time it is convenient to use basis vectors e x i and e y > shown in Fig. 4. Then, 
the six states 0, 1, 2, 3, 4 and 5 are represented respectively by 2e x > + e y / + r, 2e x / + 2e y i + r, 
e x i + 2e y i + r, e y > + r, r and e x > + r with r = e x > + e y >. The conditions in terms of 
isumx-2*isumy, 2*isumx-isumy and isumx+isumy will be used for the determination of 
the block six-clock state. 

4. Summary and Discussions 

We have presented the multi-spin coding for the Monte Carlo simulation of the three- 
state Potts model. The multi-spin coding technique has also been applied to the block-spin 
transformations of the Potts model both with the F and AF order parameters. The emphasis 
should be put on the fact that we have totally excluded if -statement in our code; as a result, 
the code can perform ultrafast computations. 

As we have mentioned in the Introduction, some earlier applications of the multi-spin 
coding for the Ising model §i'@ as well as a recent one! showed how to store several spins 
belonging to the same lattice in a word. This type of coding is effective for saving the 



9 



Fig. 4. Possible vectors obtained by the vector sum of eight six-clock spins. 



memory of computers. But not so dramatic acceleration in the simulation speed is expected. 
On the other hand, the method we have presented is based on yet another type of the 
multi-spin coding, which was originally proposed by Bhanot et al. for the Ising modela 
and modified by MichaelJHJ In these methods and the present one as well, the multi-spin 
coding technique is used in order that a number of independent systems are simulated 
simultaneously. These simultaneous simulation methods do not help reducing the memory 
consumption; but we can expect a great reduction in the computation time. Since the 
memory restriction is getting more and more relaxed these days, the simultaneous simulation 
technique is advantageous over the earlier multi-spin coding methods. A simultaneous 
simulation method using 3 bits/spin coding has also been used for the Ising model a decade 
ago.E3 

A short comment will be made here on the performance of the present code. We have 
achieved the speed of 0.81 G Potts-spin flips per second without any measurement in the 
case of 32-bits multi-spin coding of the 3D Potts model for periodic boundary conditions 
on the HITAC S-3800/480 system at the Computer Center of University of Tokyo. This is 
comparable to the speed of 2.10 G spin flips per second for the 3D Ising model on the same 
machine: For the 3D Potts model the number of the logical commands for a spin update is 
110, which is about four times of that for the 3D Ising model. With the measurement of the 
energy and the magnetization, the performance for the 3D Potts model is changed to 0.49 
G spin flips per second: When one makes block-spin transformations three times with the 
measurement of block-spin magnetizations, the speed becomes 0.44 G spin flips per second. 
We should note that the present code can be extended to the case of four-state Potts model, 
five-state Potts model, and so on, although the coding will be more cumbersome. 

Using this fast code, we have investigated several interesting topics of the Potts model. 
The order-parameter distribution functions of the F three-state Potts model has been stud- 
ied in 2D and 3D@ where attention has been paid to the vector character of the order 
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parameter. We have also treated the 3D AF Potts model for a careful study of the critical 
phenomena, using the MCRG method combined with finite-size scaling analysis; we have 
confirmed that the 3D AF three-state Potts model is in the XY universality class within 
statistical errors 

Another problem is the random three-state Potts model in 3D. We have studied the 
± J model with general asymmetric probability weights. The overall phase diagram of the 
random Potts model as a function of the concentration of coupling and the temperature 
has been obtained. In the AF-rich region, we have paid attention to the universality of 
critical phenomena. We have also shown the rounding of the first-order transition due to 
the randomness in the F-rich region. 

The details of the physical results are given in each publication^ 
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Appendix A 

c 

c Multi-Spin coding of the 3-dimensional 3-state Potts model 
c with ferromagnetic (F) , antif erromagnetic (AF) 

c and random couplings 

c 

subr out ine sblt ce ( i spO , i spl , j x , j yr , j yl , j zr , j zl , ira , irb , iengd) 
parameter (nx=64 , nx2=nx/2 , nxy2=nx2*nx , nla2=nxy2*nx) 
parameter (nword=2**24) 

c ommon/ j xt ab/ j xt bO ( : nwor d- 1 ) , j xtb 1 ( : nword- 1 ) 
& , jxtb2(0:nword-l) 

dimension isp0(0:nla2-l,0: 1) ,ispl(0:nla2-l,0: 1) 

dimension jx(0:nla2-l) , jyr(0:nla2-l) , jyl(0 :nla2-l) 
& , jzr(0:nla2-l) , jzl(0:nla2-l) 

dimension ira(0:nla2-l) , irb(0:nla2-l) 

dimension iengd(0:nla2-l,0:3) 

do 10 la=0,nla2-l 

c old spin -> isold 

isoldO = isp0(la,0) 
isoldl = isp0(la,l) 

c set trial new spin -> isnew 

iscoin = ira(la) 

isnewO = iand(not (isoldO) , iscoin) 
isnewl = not (ior (isoldl , iscoin) ) 

c energy difference for spin 1 -> idel 
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neq : flag to tell whether n.n. spins are 

the same (0) or not (1) 
ide : e(old) - e(new) + 1 

neqold < neqnew then ide = (2) for F (AF) 

neqold = neqnew then ide = 1 

neqold > neqnew then ide = 2 (0) for F (AF) 

Note that (neqold = 0) and (neqnew = 0) are not realized 
at the same time. 

isO = ispl(la,0) 
isl = ispl(la,l) 

neqold = ior(ieor(isoldO,isO) , ieor (isoldl , isl) ) 

neqnew = ior(ieor(isnewO,isO) , ieor (isnewl , isl) ) 

ide 10 = iand (neqold, neqnew) 
for F : 

ide 11 = not (neqnew) 
for AF : 

ide 11 = not (neqold) 
for random : 

idell = ior(iand( jbond, not (neqnew) ) ,not(ior(jbond, neqold))) 

energy difference for spin 2 -> ide2 

isO = ispl(jx(la) ,0) 
isl = ispl (jx(la) , 1) 

neqold = ior(ieor(isold0,is0) , ieor (isoldl , isl) ) 
neqnew = ior(ieor(isnew0,is0) , ieor (isnewl , isl) ) 
ide20 = iand (neqold, neqnew) 
ide21 = not (neqnew) 
ide21 = not (neqold) 

ide21 = ior(iand( jbond, not (neqnew) ) ,not(ior(jbond, neqold))) 

energy difference for spin 3 -> ide3 

isO = ispl(jyr(la) ,0) 
isl = ispl (jyr (la) , 1) 

neqold = ior(ieor(isold0,is0) , ieor (isoldl , isl) ) 
neqnew = ior(ieor(isnew0,is0) , ieor (isnewl , isl) ) 
ide30 = iand (neqold, neqnew) 
ide31 = not (neqnew) 
ide31 = not (neqold) 

ide31 = ior(iand( jbond, not (neqnew) ) ,not(ior(jbond, neqold))) 
energy difference for spin 4 -> ide4 



13 



isO = ispKjyl(la) ,0) 
isl = ispl(jyl(la) ,1) 

neqold = ior (ieor (isoldO , isO) , ieor (isoldl , isl) ) 
neqnew = ior (ieor (isnewO , isO) , ieor (isnewl , isl) ) 
ide40 = iand (neqold, neqnew) 
ide41 = not (neqnew) 
ide41 = not (neqold) 

ide41 = ior (iand(jbond, not (neqnew) ) ,not(ior(jbond, neqold))) 



energy difference for spin 5 -> ide5 



isO = ispl (jzr (la) ,0) 
isl = ispl (jzr (la) , 1) 

neqold = ior (ieor (isoldO , isO) , ieor (isoldl , isl) ) 
neqnew = ior (ieor (isnewO , isO) , ieor (isnewl , isl) ) 
ide50 = iand (neqold, neqnew) 
ide51 = not (neqnew) 
ide51 = not (neqold) 

ide51 = ior (iand(jbond, not (neqnew) ) ,not(ior(jbond, neqold))) 



energy difference for spin 6 -> ide6 



isO = ispKjzl(la) ,0) 
isl = ispl(jzl(la) ,1) 

neqold = ior (ieor (isoldO , isO) , ieor (isoldl , isl) ) 
neqnew = ior (ieor (isnewO , isO) , ieor (isnewl , isl) ) 
ide60 = iand (neqold, neqnew) 
ide61 = not (neqnew) 
ide61 = not (neqold) 

ide61 = ior (iand(jbond, not (neqnew) ) ,not(ior(jbond, neqold))) 



sum of energy differences for 2 spins 



isumaO = ieor (idelO , ide20) 

isumal = ior (ieor (idell , ide21) , iand(idelO , ide20) ) 
isuma2 = iand(idell , ide21) 
isumbO = ieor (ide30 , ide40) 

isumbl = ior (ieor (ide31 , ide41) , iand(ide30 , ide40) ) 
isumb2 = iand(ide31 , ide41) 
isumcO = ieor (ide50 , ide60) 

isumcl = ior (ieor (ide51 , ide61) , iand(ide50 , ide60) ) 
isumc2 = iand(ide51 , ide61) 



isuma = isuma + 2 
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isuma2 = ior (isumal , isuma2) 
isumal = not (isumal) 

c isumw = isuma + isumb 

isumwO = ieor (isumaO , isumbO) 
icr = iand(isumaO, isumbO) 
iwk= ieor (isumal , isumbl) 
isumwl = ieor (iwk, icr) 

icr = ior (iand(isumal , isumbl) , iand(iwk, icr) ) 
iwk= ieor(isuma2,isumb2) 
isumw2 = ieor (iwk, icr) 

isumw3 = ior (iand(isuma2 , isumb2) , iand(iwk, icr) ) 
c isum = isumw + isumc 

c isum : total trial energy difference (old-new) +8 
c (2 <= isum <= 14) 

isumO = ieor(isumwO,isumcO) 
icr = iand(isumwO, isumcO) 
iwk= ieor (isumwl , isumc 1) 
isuml = ieor (iwk, icr) 

icr = ior (iand(isumwl , isumcl) , iand(iwk, icr) ) 
iwk= ieor(isumw2,isumc2) 
isum2 = ieor (iwk, icr) 

icr = ior (iand(isumw2, isumc2) , iand(iwk, icr) ) 
isum3 = ior (isumw3, icr) 

c jxtb : X-table for Boltzmann factor 

c (0 <= jxtb <= 6) 

c iflip : flag for spin flip 

c isum + jxtb >= 8 then flip 

ir = irb(la) 
jtO = jxtbO(ir) 
jtl = jxtbl(ir) 
jt2 = jxtb2(ir) 

icr = iand(ieor (isuml , jtl) ,iand(isumO, jtO)) 
icr = ior (iand(isuml , jtl) , icr) 
icr = iand(ieor(isum2, jt2) ,icr) 
icr = ior(iand(isum2, jt2) ,icr) 
iflip = ior (isum3, icr) 

c new spin -> ispO(la) 
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ispO(la,0) = ior (iand(if lip, isnewO) , iand(not (if lip) , isoldO) ) 
ispO(la,l) = ior (iand(if lip, isnewl) , iand(not (if lip) , isoldl) ) 

c energy difference (old-new) +8 (2 <= iengd <= 14) 

iengd(la,0) = iand(if lip, isumO) 

iengd(la,l) = iand(if lip, isuml) 

iengd(la,2) = iand(if lip, isum2) 

iengd(la,3) = ior (not (if lip) , isum3) 

10 continue 

end 
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